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We present an efficient new computational method for calculating the binding energies of the 
bound states of ultracold alkali-metal dimers in the presence of magnetic fields. The method is based 
on propagation of coupled differential equations and does not use a basis set for the interatomic 
distance coordinate. It is much more efficient than the previous method based on a radial basis set 
and allows many more spin channels to be included. This is particularly important in the vicinity 
of avoided crossings between bound states. We characterize a number of different avoided crossings 
in Cs2 and compare our converged calculations with experimental results. Small but significant 
discrepancies are observed in both crossing strengths and level positions, especially for levels with I 
symmetry (rotational angular momentum L = 8). The discrepancies should allow the development 
of improved potential models in the future. 



I. INTRODUCTION 

Ultracold Cs atoms are of great interest for a num- 
ber of experiments, which haveproduced a Bose-Einstein 
condensate of such atoms formed a cold cloud 
of Cs2 dimer molecules @], probed three-body Efimov 
physics studied collisional shifts Q or quantum 
scattering [f| of atomic clock states, carried out high- 
resolution molecular spectroscopy [|| or used magnetic 
fields to switch among a variety of very weakly bound 
molecular states of the Cs2 dimer 0, [1[ . These experi- 
ments all depend upon and take advantage of the colli- 
sional interactions between two Cs atoms. Consequently, 
accurate theoretical and computational models of near- 
threshold Cs atom scattering and bound states are neces- 
sary for maximum understanding of existing experiments 
and for making quantitative predictions for new experi- 
mental domains. 

Because of the complex spin structure of two ground- 
state Cs atoms, many different near-threshold bound 
states exist and have different magnetic moments. They 
thus tune differently with magnetic field. When one of 
these bound states crosses a collision threshold, a low- 
energy scattering resonance occurs, commonly known 
as a Feshbach resonance. Extensive study of such res- 
onances has allowed the construction of quite accu- 
rate coupled-channel models for calculating the magnetic 
field-dependent scattering and bound-state properties 
near collision thresholds @, H, 0, E3, EH E2 ■ These mod- 
els incorporate the electron and nuclear spins, their mu- 
tual interactions, and the adiabatic Born-Oppenheimer 
potentials for the X 1 E+ and a 3 E+ molecular states that 
correlate with two 2 S 1( /2 ground-state Cs atoms. By ad- 
justing the model parameters to fit the measured mag- 
netic fields for resonances in different scattering chan- 
nels, the model quite accurately predicts near-threshold 
scattering properties and the binding energies of weakly 



bound states within a few GHz of threshold. Such thresh- 
old models can also be adapted to treat three-body inter- 
actions, for which an accurate knowledge of the threshold 
two-body bound states is necessary [13| . The models are 
sensitive to relatively few parameters, and may or may 
not be adequate when extended into new experimental 
domains. 

Recently, Mark et al. [£, 8] have characterized a num- 
ber of avoided crossings between levels bound by only 
E/h~ 5 MHz with respect to the energy of two separated 
Cs atoms in their lowest-energy Zeeman sublevels. Using 
time-dependent magnetic field ramping, they were able to 
convert two Cs atoms into a number of different molec- 
ular states with different rotational quantum numbers 
and magnetic moments. Most of the bound states are 
well described by the existing coupled-channel model in 
regions far from avoided crossings. However, character- 
izing the avoided crossings themselves presents problems 
for the existing computational methods. In particular, 
Ref. [T2| calculated bound states using a method based 
on a basis set expansion of the radial wavefunctions in 
a discrete variable representation (DVR). This method 
can use only a restricted spin basis in determining the 
molecular bound states because of the large number of 
grid points required. 

The present paper develops an improved computa- 
tional method that is necessary to calculate and un- 
derstand the avoided crossin gs i n Cs2- This method 
uses a propagator approach [14| in place of a radial 
basis set to represent the molecular bound states. It 
can readily be adapted to threshold states of other 
molecules [lot fl(| . The propagator approach is compu- 
tationally much cheaper than the DVR approach and as 
a result can include many more coupled spin channels. 
The new approach is used to compare the calculated and 
observed properties of the avoided crossings, in order 
to identify aspects of the ground-state coupled-channel 



model for CS2 that are still in need of improvement. 



II. COMPUTATIONAL METHODS 

The present work solves the bound-state Schrodingcr 
equation for CS2 using two independent methods. In ci- 
ther case the Hamiltonian may be written 
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+ hi + h 2 + V(R), (1) 



where /i is the reduced mass and L 2 is the operator for 
the end-over-end angular momentum of the two atoms 
about one another. The monomer Hamiltonians includ- 
ing Zeeman terms are 



hj = (ij 



j- 



(2) 



where si and s 2 represent the electron spins of the two 
atoms and i\ and 12 represent nuclear spins. g e and g n 
are the electron and nuclear <?-factors, is the Bohr 
magneton, and s z and i z represent the z-components of 
s and 1 along a space-fixed Z axis whose direction is de- 
fined by the external magnetic field B. The interaction 
between the two atoms V(R) is given by Stoof et al. [Tt} 
as the sum of two terms, 



V(R) = V C (R) + V d (R) . 



(3) 



Here V°(R) = V (R)V {0) +V 1 (R)V^ 1) is an isotropic po- 
tential operator that depends on the potential energy 
curves Vq(R) and V\{R) for the respective X lv J+ singlet 
and a 3 E+ triplet states of the diatomic molecule. The 
singlet and triplet projectors and project onto 
subspaccs with total electron spin quantum numbers 
and 1 respectively Figure [T] shows the two potential en- 
ergy curves for CS2. The V d (R) term represents small, 
anisotropic spin-dependent couplings that are responsi- 
ble for the avoided crossings discussed in this paper and 
are discussed further in Section HTT1 below. 

The first method for finding eigenvalues is a conven- 
tional full matrix diagonalization in a discrete variable 
representation (DVR) [lH. It uses a basis set made up 
of products of internal and radial functions. The inter- 
nal Bosc-symmctrized basis set is made up of functions 
in which the operators L 2 and hj are diagonal, that is, 



\a 1 mf 1 )\a 2 mf2)\LM L }, 



(4) 



where \LMl) and \ajirifj} respectively represent the 
eigenstates of L 2 and the B-dependent monomer Hamil- 
tonian hj, and where Ml and rrifj are projection quan- 
tum numbers along the magnetic field direction. When 
i? = 0, \ajirifj) = \(sjij)fjtnfj), where fj is the total 
spin of atom j and mfj is its space- fixed projection. As 
B increases from zero, different fj values become mixed. 
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FIG. 1: Molecular potential energy curves Vo (R) and Vi(R) 
for the respective singlet and triplet states of CS2 correlat- 
ing with two separated 2 Si/2 ground-state atoms. The inset 
shows an expanded view of the long-range potentials separat- 
ing to the two different / = 3 and 4 hyperfine states of the 
2 Si/2 atom with nuclear spin i=7/2 and magnetic field B — 0. 
The inset shows the adiabatic potentials obtained from diag- 
onalizing the matrix form of the operator hi + hi + V(R) 
at each R for the case of L = 0, Mf = +6. There are 5 
channels, and the 3 + 4 and 4 + 4 separated-atom limits are 
doubly degenerate at B = 0. All 5 channels have the same 
long-range variation as —Ce/R 6 , with Ce = 6860 Eha® [l2l | 
{E h = 4.3597X 10~ 18 J is the Hartree and a =0. 0529177 nm is 
the Bohr radius). The level crossings discussed in this paper 
are for very weakly bound levels that lie within about E/h ~ 5 
MHz of the dissociation limit to two {/w/} = {3, +3} atoms 
in the magnetic field range from mT to 5 mT. 



The DVR radial functions are unevenly spaced colloca- 
tion points obtained from a nonlinear coordinate trans- 
formation (l9j . 

This DVR method requires diagonalizing a large NxN 
matrix, the dimension of which is given by the product 
of the number of spatial collocation points jV c and the 
number of spin basis functions N s . We use the LAPACK 
subroutine DSPEVX to find a selected range of eigen- 
values and eigenvectors [2(|. In order to use a direct 
diagonalization procedure to calculate the bound-state 
energies [2l| shown in Rcfs. [1, [HI, the magnitude of 
N = N C N S was limited to around 25000 using a proces- 
sor with 4 GB of memory. With N c « 800, in order to 
give 5 points per node with about 150 nodes for threshold 
wave functions, the number of spin basis functions is thus 
restricted to be about N s w 35. When this is fewer than 
is needed for a complete calculation, an approximation 
scheme becomes necessary, as described in Section [IHI 

The second method avoids the use of a basis set for the 
interatomic distance R and instead relies on propagation 
of coupled differential equations In this case the 

Bose symmetrized basis set used is a fully decoupled set, 
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The compound channel index k is used to simplify nota- 
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tion and implies values of all the quantum numbers in the 
basis set. While the choice of the basis sets in Eqs. ^ 
and ([5|) represent different approaches, they are equiv- 
alent for representing molecular energy levels when the 
two basis sets span the same space. There is a simple 
unitary transformation between the two basis sets. The 
matrix elements of the different terms in the Hamiltonian 
in basis set ((SJ) are given in the Appendix. 

In the propagation method, we expand the total wave- 
function for state n as 



(0) 



Substituting into the Schrodingcr equation and project- 
ing onto each channel function in turn gives a set of cou- 
pled equations for the radial channel functions tpkn{R), 
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Y,W jk {R)-e5 jk ]^ kn {R), (7) 



where 6j k is the Kronecker delta, e is the energy E scaled 
by 2/j,/h 2 , and 
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h 2 + V(R) 



dr, 



(8) 

where dr indicates integration over all coordinates except 
R. If there are N s basis functions, the required solu- 
tion ip n (R) is a column vector of order N s with elements 
ipkn{R)- However, Eq.[7]has N s independent solution vec- 
tors at any energy, so that until the boundary conditions 
are applied tp n (R) is an iV s x N s wavefunction matrix. 

The Schrodinger equation can be solved to find an 
N s x N s wavefunction matrix at any energy E. In practice 
it is numerically stabler to propagate the log-derivative 
matrix Y(R) = [dip n / 'di?][^ n (i?)] _1 . However, a solu- 
tion that satisfies bound-state boundary conditions can 
be found only at the eigenvalues E n . Solutions are propa- 
gated outwards from a point iJ m j n in the inner classically 
forbidden region and inwards from a point i? max at long 
range to a matching point i? m id- The outwards and in- 
wards solutions are designated Y + (R) and Y~ (R). If E is 
an eigenvalue of the coupled equations, there must exist a 
wavefunction vector ip n (R mid ) = ^(-Rmid) = ^(Rmid) 
for which the derivatives also match, 



dR 



dip; 



dR 



(9) 



so that 



Y + (R mid )MRmid) = Y~{R mid )MRmi<i)- (10) 

Thus ipn(Rmid) is an eigenvector of Y + (R mid )-Y~ (R mid ) 
with eigenvalue 0. It is thus possible to locate eigen- 
values of the Schrodingcr equation by propagating solu- 
tions of the coupled equations and searching for zeroes 
in the eigenvalues of the log-derivative matching matrix 



Y + (R mid ) — Y~(R m i d ) as a function of energy. This ap- 
proach is much stabler for large multichannel problems 
than the older approach [22j of searching for zeroes of the 
determinant of the matching matrix. 

The major advantage of the propagator method is that 
the matrices handled are only of dimension iV s x 2V 8 , 
where N s is the number of internal basis functions. The 
computational cost is proportional to -/V s 3 but only linear 
in the number of propagation steps. By contrast, a full 
diagonalization with N c radial basis functions (colloca- 
tion points) involves matrices of dimension N S N C x N S N C . 
The computational cost is proportional to N^N^. Since 
N c typically needs to be greater than 500 for the present 
application, the propagator approach is much cheaper. 

The BOUND program [23| is a general-purpose pack- 
age to solve the bound-state Schrodinger equation using 
propagator methods. The algorithms used are described 
in more detail in Ref. [lj|. For the purpose of the present 
work we have generalised the BOUND package in three 
significant respects: 

1. We have generalised the structure of the code so 
that it can handle coupled equations in basis sets 
that are not diagonal at R — oo; 

2. We have implemented the specific set of coupled 
equations required for Cs2 with the basis set of Eq. 

El 

3. We have added an option to use the log-derivative 
propagator of Alexander and Manolopoulos [24| . 
which is based on Airy functions and allows very 
large step sizes at long range. 

In the presence of a magnetic field, the only rigorously 
conserved quantum numbers are M to t = mfi + m fi + 
Mi = m s i + m S 2 + m%i + ma, + Ml and the total par- 
ity (— 1) L . This leads to an infinite number of channels. 
However, L and Mp = m/i + m/2 are very good ap- 
proximate quantum numbers because the only term in 
the Hamiltonian that is off-diagonal in them is the small 
anisotropic coupling term V A . In either computational 
approach it is possible to restrict the number of chan- 
nels by selecting only one or a few values of L and all 
or a subset of possible Mp values. Here we consider 
the case studied experimentally by Mark et al. [§}, who 
used Cs atoms is their lowest energy hyperfine state with 
rrif = +3 to make Cs2 molecules with M tot = +6. The 
number of channels with L = 0, 2, 4, 6, 8, including all 
allowed Mp values, arc 5, 23, 46, 76, 103, respectively. 
Thus, for example, a full calculation including all chan- 
nels with L = 4, 6 and 8 requires 225 channels. 

In practical terms, for example, a run with the DVR 
method to find 28 bound states within 3 GHz of the 
E = threshold for Cs2 for a single magnetic field with 
30 channels and 720 collocation points took about 7 hours 
on an 2.4 GHz processor. With the propagator approach 
we were able to find selected near-dissociation levels for 
30 channels in about 40 seconds per level with a 2.0 GHz 



processor. The great advantage of the propagator ap- 
proach was demonstrated by our ability to find levels 
with 225 channels in about 45 minutes per level. A cal- 
culation with 225 channels would not be possible at all 
using the DVR method with a direct eigenvalue solver. 



III. COMPARISON OF COMPUTATIONAL 
RESULTS 

The DVR and propagator calculations described here 
both use the same potential model, with the parameters 
given by Chin et al. [l2j . The potential energy curves are 
based on the ab initio calculations of Krauss and Stevens 
[25| . The singlet and triplet scattering lengths as and 
(Xt, the long-range coefficients Cq and C%, and a scaling 
factor Sc for the second-order spin-orbit coupling were 
adjusted by Chin et al. to reproduce a substantial number 
of Feshbach resonances with L < 4. 

Figure [5] shows an example of weakly bound levels of 
the CS2 molecule with M to t = +6 in the mT to 6 
mT range of B. Many of these levels have been probed 
in the experiment of Mark et al. Q. The figure also 
shows the bound-state classification scheme of Chin et 
al. [Hj], namely FL(Mp), where F is the resultant of the 
separated-atom spins f\ and f% and Mp is its projection 
defined above. Like fi and /2, F is a good approximate 
quantum number for labeling near-threshold levels at low 
B. Quantum numbers L = 0, 2, 4, 6, 8 are represented by 
labels s, d, g, i, I, respectively. Ml need not be specified 
since Ml = M tot — M F . Fig. [5] shows levels with L < 4 
obtained from a DVR calculation that included only ba- 
sis functions for a single L and Mp. This neglects the 
small off-diagonal couplings between levels with different 
L and Mp quantum numbers due to V d , so that levels 
of different symmetry show crossings rather than avoided 
crossings in Fig. [2] 

For ground-state alkali-metal atom interactions, the 
V d operator has the form of spin-dipolar coupling 

V d (R) = X(R) (s x • h - 3(s x • e R )(s 2 ■ e R )) , (11) 

where €r is a unit vector along the internuclcar axis and 
A is an i?-dependent coupling constant, which for our 
model is 



6g(2) 



\{R) = EW 



1 



(R/ao) 3 



0.071968e 



-0.83[(.R/a )-10] 



(12) 

where a w 1/137 is the fine structure constant. At large 
R the coupling becomes the long-range dipolar interac- 
tion between the spins on the separated atoms that varies 
as l/i? 3 [Tt], H(| . In the short-range region of chemical 
bonding the magnitude of X(R) is primarily determined 
by the second-order spin-orbit coupling term represented 
by the exponential term [l^, H3, Ha [29| . 

The crossings in Fig. [2] become avoided crossings when 
the small interactions due to V d are taken into account. 
The energy splitting at the crossing varies greatly, de- 
pending on the quantum numbers of the two levels. In 
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FIG. 2: Bound state energy E/h as a function of B for levels 
of the Cs2 molecules with even L < 4 and M to t = +6. Ener- 
gies are given relative to the energy of two Cs atoms in their 
ground Zeeman sublevel (/ = 3, mf = +3). The FL(Mf) la- 
beling scheme is shown for each level. Off-diagonal coupling 
between levels with different FL(Ml) quantum numbers is 
neglected in this calculation. 



first order, the V d operator couples states FL(Mp) and 
F'L'(M' F ) according to the selection rules \L — L'\ = or 
2, \F-F'\ = 0, 1, or 2, and \M F -M' F \ = 0, 1, or 2. These 
selection rules immediately follow from the tensor form of 
the operator in Eq. (fTTj) , as given by Stoof et al. [lj} , who 
write Eq. (jlip as a sum of products of L = 2 spherical 
harmonic components Ylm l {gr) and rank 2 spin tensor 
components. We refer to a crossing as direct when there 
is a first order coupling of the two states involved and 
indirect when there is not. 

The success of a calculation of the CS2 energy levels 
and their avoided crossings depends on the sufficiency of 
the basis set expansion of the wave function. Suppose we 
wish to calculate the energy of one FL(Mp) state that 
crosses a different F'L'{M' F ) state. It is necessary to 
include sufficient basis functions to represent each state 
adequately, and to represent their interaction. This is 
simplified by taking advantage of the selection rules de- 
scribed above. In order to represent a level with a given 
FL(Mp), it is necessary to include all basis functions 
with the same set of three quantum numbers, since such 
levels are coupled by terms due to the strong central po- 
tential V c . A level calculated with such a basis is coupled 
through the V d operator to other levels in which one or 
more of the three quantum numbers are different. Such 
off-diagonal coupling causes shifts in level positions and 
also induces avoided crossings. 

In the propagator calculations, the basis set usually in- 
cludes all functions with L and V of the levels in question 
consistent with M tot . Additional basis functions with 
different quantum numbers Li are added to account for 
shifts and crossings due to coupling of L or L' with Li. 
The propagator basis is specified by giving L, L' , and a 
list of additional values Li needed to account for higher- 



order coupling. In the DVR calculations, the basis sets 
are additionally limited by restricting the calculation to 
functions with L(M F ), L'(M' F ) and additional quantum 
numbers Li(MFi) as needed. Thus the basis set is spec- 
ified by giving the list L(M F )L'(M F )[Li(M F ,i)]. Some 
propagator calculations were done with a similarly re- 
stricted list to verify that the two methods gave exactly 
equivalent results. Neither the propagator nor DVR cal- 
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FIG. 3: Example of coupling between different L(M F ) sym- 
metry blocks with the symmetry of the dipole-dipole inter- 
action of V d . Each block represents a Hamiltonian matrix 
for spin states with the L(Mf) values indicated. The labels 
"x" and "y" indicate the existence of nonvanishing coupling 
due to V ; a "0" indicates no coupling. The case shown is 
for a g(&) and a g(S) level, which have no direct coupling. 
The left panel shows the symmetries that give rise to second- 
order interactions between the two levels and thus contribute 
to the strength of the avoided crossing between them. The 
right panel shows a truncated set of interactions through in- 
termediate d(4) and d(5) levels. 



Figure [3] illustrates the size of the basis set needed, as 
governed by the selection rules on V coupling. Since 
the matrix elements of V are relatively small, they are 
normally of practical significance only through second 
order. Thus it is necessary to include only intermediate 
levels with Lj and M F ,i that differ from L or L 1 and M F 
or M' F by at most 2 units. Any higher-order couplings 
would be much smaller than those discussed here. Thus, 
in order to represent the crossing of a 6g(6) and a 4^(3) 
level, for which there is no first-order direct coupling, d-, 
g-, or i-basis functions with M F — 4 and 5 need to be 
included in the basis, as shown in Fig. [3] To represent 
additional second-order shifts of the two g levels, d-, g- 
and i-basis functions with M F = 1, 2, 7 and 8 also need 
to be added. 

Figure [4] illustrates calculations with different basis 
sets, comparing energies calculated with the propagator 
and DVR methods for the crossing of the 4<?(3) and 65(6) 
levels near 1.0 mT. Table U tabulates the positions and 
strengths of this crossing, as well as a number of oth- 
ers. The position Bq of the crossing is defined as the 
field at which the two levels are closest together and the 
strength 2V is the minimum of the difference between 
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FIG. 4: Calculated energy levels E/h with M tot = +6 as a 
function of magnetic field B near the crossing of the 4g(3) 
level with the 6g(6) level near 1.0 mT. The points and solid 
line show the propagator calculation with a sdgi basis set. 
The dashed lines show the crossing levels from two uncoupled 
DVR calculations with g(3) or g(Q) basis functions only. The 
dash-dot and dotted lines show the crossing levels from DVR 
calculations with added g(4, 5) and rf(4, 5) functions respec- 
tively. The DVR calculation with i(4, 5) basis functions is not 
shown, but lies near the uncoupled crossing and has a very 
small splitting, indicating very weak second-order coupling 
through distant i states. 



the two energies as a function of B; 2V is used since 
the splitting is twice the effective coupling matrix ele- 
ment V in a 2-level representation of the crossing 
We have verified that the two methods give identical re- 
sults within numerical accuracy when exactly equivalent 
basis sets are used. Since there is no direct interaction 
between the two crossing levels in this case, the splitting 
at the crossing originates principally in second-order in- 
teractions mediated through distant levels of d, g, or i 
symmetry with M F = 4 or 5. However, as mentioned 
above, second-order couplings to levels with other M F 
values can cause additional shifts. Both bound and scat- 
tering states can contribute, and the contribution from 
any given distant state varies inversely with the its sepa- 
ration in energy from the crossing. Intermediate g levels 
are the closest in energy to the crossing, whereas interme- 
diate i levels are the most distant. In Figure 31 the sdgi 
basis set used in the propagator calculation is effectively 
complete. It may be seen that a calculation including 
only the c/(4, 5) intermediate states captures most of the 
crossing strength but does not reproduce the level shifts 
well. Conversely, a calculation including only the d(4, 5) 
states gives a crossing strength that is much too small 
but overestimates the level shifts. The contributions to 
the crossing strength from different intermediate states 
are far from additive. There are no experimental results 
for this crossing. 

Figure [5] illustrates a different case, a 4g(4) — 65(6) 
crossing with a splitting that is about 8 times larger at 
the crossing. This is a case where the two states in- 
volved have a direct coupling to one another through 




FIG. 5: Calculated energy levels E/h with M to t = +6 as a 
function of magnetic field B near the crossing of the 4<?(4) 
level with the 6g(6) level near 1.0 mT. The points and solid 
line show the propagator calculation with a dgi basis set. The 
dashed lines shows the crossing levels from two uncoupled 
DVR calculations with gr(4) or g(6) basis functions only. The 
dash-dot lines show the avoided crossing from a DVR calcu- 
lation with only the direct coupling in the g(4, 6) basis set 
included. The doubled-headed arrow at the position of the 
propagator crossing shows the measured splitting The 
actual experimental crossing was observed 0.046 mT lower in 
B value than the propagator crossing. 



V d . While additional second-order coupling can change 
the position and strength of the crossing slightly, the di- 
rect coupling is dominant and the restricted-basis DVR 
calculation agrees much better with the full propagator 
calculations. Both are in reasonable agreement with the 
measured splitting of the crossing However, the cal- 
culated position in B needs to be shifted by —0.046 mT 
to agree with the measured position Q. This remaining 
discrepancy reflects a real deficiency in the parameters of 
our potential model as discussed below. 

Figure [5] shows the difference between the upper and 
lower branches of the crossing for the case of the 4<f (4) — 
65(6) crossing near 4.5 mT. This is another case of direct 
coupling, where the measured and calculated crossings 
agree well in coupling strength, although the calculated 
position needs to be shifted by +0.034 mT to agree with 
the measured one. 

Table U show comparisons between the propagator 
and DVR calculations for a number of other crossings 
in Figure [5] Crossing positions generally agree within 
about 0.01 mT among the different basis sets. Relatively 
good agreement between propagator and DVR coupling 
strengths 2V is seen in the cases where there is direct cou- 
pling between the two crossing levels, or where the DVR 
method includes all second-order intermediate states al- 
lowed by the symmetry of the V d operator. However, 
for higher- L crossings it was usually necessary to select a 
subset of the allowed intermediate states to make DVR 
calculations feasible. In such cases the DVR method can 
give unreliable results, depending on the choice of rc- 



FIG. 6: Calculated energy difference A/h between the 6g(6) 
and 4d(4) levels with Mtot = +6 as a function of magnetic 
field B. The solid line is from a propagator calculation with 
the sdgi basis. The diamonds show the experimental results 
obtained by Ferlaino el al. using their more accurate field 
modulation method [3(J. The dashed line shows the calcu- 
lated points shifted by +0.034 mT. The DVR calculation (not 
shown) with direct coupling included in the d(4)g(6) basis is 
virtually identical to the dashed line when the DVR results 
are shifted by +.062 mT. 



stricted basis set. 

Figure [7] shows calculated bound states for s and d 
levels on a broader energy scale. (Levels with other sym- 
metries are not shown). A DVR calculation with a full 
sd basis is possible in this case. The 6s(6) and 4d(4) 
uncoupled levels show two crossings. The low-field cross- 
ing around 0.24 mT occurs near the observed location 
(0.72 mT) of a three-body Borromean state of the CS3 
trimer associated with the exotic Efimov physics of this 
species 0]. Lee et al. [H| used the last two 6s(6) two- 
body states of the CS2 dimer to construct the parame- 
ters for full three-body calculations of bound states and 
recombination coefficients in the mT to 3 mT range. 
While their method was able to give semi-quantitative 
agreement with the measurements, the avoided crossing 
of the 6s(6) level with the 4<f(4) level needs to be taken 
into account in subsequent calculations because of the 
mixed spin character of the target molecular state pro- 
duced by the three-body recombination in this region of 
B. The strong s — d interactions modify the s-wave scat- 
tering length at small B, but this is easy to take into 
account by including s and d basis functions in scatter- 
ing calculations. 

The higher-field 6s(6) — 4d(4) crossing near 4.8 mT 
has been studied in Rcfs. [I,l3l[. Fi gurc [9] shows an 
expanded view of the very-near-thrcshold region of this 
crossing and the additional 6s(6) — 2<?(2) crossing near 
5.4 mT. The interaction between s and d states results 
in an overall shift in the binding energy of the 6s(6) level, 
where the uncoupled level is too high in energy. This case 
illustrates one advantage of the propagator method over 
the DVR method. The latter has to use a finite range of 
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spatial points and is restricted by the length of the "box" 
in which the calculation is carried out. When this length 
is too large, the number of the spatial collocation points 
can become too large for practical calculations. A 5000 
ciq "box" is sufficient for levels with binding energies on 
the order of 40 kHz, since the scattering length, which 
gives an indication of the "size" of the weakly bound 
molecular levels [l5|, is on the order of 1000 cio <C 5000 
ao- Such restrictions on spatial grid do not apply to the 
propagator method, which is capable of calculating levels 
arbitrarily close to E = 0, as long as the propagation is 
to sufficiently large distances. Since the propagator used 
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FIG. 7: Energy levels E/h as a function of B for the Cs2 
molecule for L — and 2 only with Mtot = +6 (g and I levels 
are not shown). The dashed lines show the DVR levels cal- 
culated with the uncoupled s(6) and d(Mp) basis sets, where 
Mf = 4, 5, or 6. The diamonds show the results of Mark et 
al. ||. The 4d(4) level crosses the 6s(6) level twice, near 0.24 
mT and 4.77 mT. The closed circles and dotted lines show the 
levels obtained from a DVR calculation with an sd basis. A 
propagator calculation with a full sdg basis shows negligible 
differences for this case. 



IV. COMPARISON WITH EXPERIMENT 

When the basis set is sufficiently large, there is good 
overall agreement between our calculations and the ex- 
perimental measurements, as already noted in relation to 
Figs. [S] and [51 Table HTl lists other examples, including the 
crossings of the 6g{6) level with the 6l(Mp) levels shown 
in Fig. [5) Since the potentials and second-order spin- 
orbit coupling of the model were originally adjusted to 
reproduce Fcshbach resonances due to zero-energy bound 
states of d and g symmetry, the positions of crossings be- 
tween s, d, and g levels tend to be accurate to within the 
model uncertainties, which are on the order of 0.05 mT 
or less [HI IH| . On the other hand, the levels of I symme- 
try, corresponding to L = 8, arc off by up to 0.5 mT, a 
much larger amount. One plausible reason for this has to 



FIG. 8: Energy levels E/h as a function of B for the Cs2 
molecule for L — 8 only with M to t = +6. The solid lines 
show the DVR levels calculated with the uncoupled 1(Mf) 
basis sets, where Mf = 0, . . . , 6. The dashed line shows the 
6c/ (6) level for which avoided crossings have been calculated 
( □ □ 
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FIG. 9: Expanded view of the crossing in Fig. [7] of the 4d(4) 
and 6s(6) levels near 4.8 mT. The long dashed line shows 
the uncoupled calculation with the s(6) and d(4) basis sets. 
The solid lines show the propagator calculations with an sdg 
basis. The upper crossing near 5.4 mT is due to a 2g(2) level. 
The open circles show DVR calculations with a full sd basis 
in a finite box of 5000 ao- The diamonds show experimental 
results of Lange et al. [32| . 



do with the large rotational energy of the 61 (Mp) levels 
that cross the 6g(6) level. The 65(6) level has the vi- 
brational character of the second 6s (6) vibrational level 
below the lowest scparated-atom limit, with about 110 
MHz of I = 4 rotational energy added. The crossing 
6l(Mp) levels, by contrast, have the vibrational charac- 
ter of the third vibrational level below the limit, with 
about 740 MHz of rotational energy added to bring them 
near threshold. More deeply bound levels with more ro- 
tational energy can have larger errors due to deficiencies 
in the model potentials. An error of only a few parts per 
1000 in the rotational energy can lead to a 0.5 mT error 
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TABLE I: Levels crossing the 6(/(6) level of the Cs2 molecule 
with a binding energy near —5 MHz relative to the energy of 
two atoms in their lowest energy hyperfine state at each B. 
The columns label the symmetry of the crossing state, the 
computational method (propagator or DVR), the L functions 
in the basis set used for the calculation (only intermediate ba- 
sis states are listed for the DVR method, since basis states for 
the two crossing states are automatically included), and the 
position Bq, energy E/h, and splitting 2V/h of each crossing. 



State Method 


Basis 


Bo 


E/h 


2V/h 






[mT] 


[MHz] 


[kHz] 


4g(2) propagator 


a 


0.7615 


-5.180 


3.6 


propagator 


dg 


0.7670 


-5.192 


1.0 


propagator 


sdgi 


0.7617 


-5.152 


1.5 


DVR 


ff(2,6)[d(4)] 


0.7745 


-5.152 


2.3 


61(3) propagator 


gil 


0.9181 


-5.135 


21.2 


DVR g(6)Z(3)[i(4,5)] 0.9112 


-5.138 


21.7 


4g(3) propagator 


sdg 


1.0105 


-5.169 


37.3 


propagator 


sdgi 


1.0024 


-5.138 


33.5 


DVR ff (3,6)[d(4,5)] 


1.0097 


-5.134 


12.1 


DVR <?(3,6)[ S (4,5)] 


0.9935 


-5.134 


32.9 


DVR 


ff(3,6)[i(4, 5)] 


0.9937 


-5.1314 


1.1 


61(4) propagator 


gil 


1.2715 


-5.088 


43.6 


4g(4) propagator 


9 


1.368 


-5.14 


264. 


propagator 


dgi 


1.375 


-5.11 


277. 


DVR 


.9(4,6) 


1.367 


-5.10 


265. 


6s(6) propagator 


sdg 


1.8648 


-5.114 


44.8 


propagator 


sdgi 


1.8664 


-5.079 


44.8 


61(5) propagator 


gil 


2.0089 


-5.056 


66.6 


61(6) propagator 


gil 


4.3211 


-4.890 


77.6 


4d(4) propagator 


dg 


4.4403 


-4.937 


55.5 


propagator 


sdgi 


4.4766 


-4.885 


57.6 


DVR 


d(4) fl (6) 


4.4485 


-4.905 


53.0 


2g(2) propagator 


a 


5.1906 


-4.887 


< 0.1 


propagator 


dgi 


5.1176 


-4.842 


9.5 



in the crossing positions for 6l(Mp) levels. 

Additional information is contained in the coupling 
strengths that govern the closest approach 2V between 
levels at avoided crossings. In the calculations, this quan- 
tity is determined largely by the second-order spin-orbit 
contribution to V d . This is a relatively poorly deter- 
mined parameter in our model and is uncertain to about 
15% 0. 

Several different experimental methods have been used 
to determine coupling strengths. Mark et al. Q used a 
method based on Stuckclbcrg intcrferomctry, which gives 
precise measurements of the energy difference between 
the two states. Mark et al. Q used a different method 
based on integrating magnetic moment values. This gives 
absolute energies for the two states (rather than just the 
difference between them) but is now believed to overesti- 
mate the coupling strengths in some cases [33| , especially 



for crossings between states with very different magnetic 
moments. Some crossing strengths were also estimated 
from a Landau- Zcner approach. Lastly, Fcrlaino et al. 
[30 ] have used a method in which transitions are induced 
by modulating the magnetic field [34|. This is the most 
precise of the different methods. 

TABLE II: Comparison of results from the best propagator 
calculation with the experimental results for selected level 
crossings with the 6g(6) state. The columns label the sym- 
metry of the crossing state, the origin of the value, the L 
functions in the basis set used for the calculation, and the 
position Bo, energy E/h, and the energy splitting 2V/h for 
each crossing. The lines labeled "Exp" show the experimental 
values. 

State method basis Bp [mT] E/h [MHz] 2V/h [kHz] 



61(3) 


propagator 


gil 


0.9181 


-5.135 


21.2 




Exp a 




1.122(2) 




32(6) 




Exp 6 




1.1339(1) 




28(2) 


61(4) 


propagator 


gil 


1.2715 


-5.088 


43.6 




Exp a 




1.550(3) 




128(26) 


4g(4) 


propagator 


dgi 


1.375 


-5.11 


277. 




Exp a 




1.329(4) 




328(60) 




Exp c 




1.357(1) 




291.4(8) 


6s(6) 


propagator 


sdgi 


1.8664 


-5.079 


44.8 




Exp c 




1.8651(3) 




58(17) 


61(5) 


propagator 


gil 


2.0089 


-5.056 


66.6 




Exp a 




2.53(1) 




126(44) 


4d(4) 


propagator 


sdgi 


4.4766 


-4.885 


57.6 




Exp a 




4.515(4) 




240(42) 




Exp c 




4.5106(3) 




78(9) 






a. 


Reference 










1). 


Reference 










c. ] 


Reference | 


30]. 





The crossing strengths for various different levels cross- 
ing the 6<?(6) level near 5 MHz are compared with the 
available experimental values in Table [II] The most reli- 
able experimental results are those from Stiickelberg os- 
cillations and magnetic field modulation [33[ for the 
6Z(3), 4g(4), 6s(6) and 4d(4) levels. The 6/(3) and 6s(6) 
levels are indirectly coupled to 65(6), and for both these 
the calculated crossing strength is about 25% lower than 
the best experimental value. The 4g(4) and 4d(4) levels 
are directly coupled to 6g(6); for the 4g(4) level the cal- 
culated crossing strength is about 5% lower than exper- 
iment, while for the 6s (6) level the discrepancy is larger 
but is within the experimental error bars. This suggests 
that the strength of the coupling term V d (R) is under- 
estimated but within the error range of Leo et al. (lOj . 

Some of the other crossings in Table HI1 show larger dif- 
ferences between experiment and theory, but in all these 
cases the experimental value was obtained using the less 
reliable magnetic moment method. The possible exper- 
imental errors for the magnetic moment approach are 
illustrated by the 4d(4) crossing, where it gives a cross- 
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ing strength a factor of 3 larger than the more accurate 
magnetic field modulation method. It would be very in- 
teresting to remeasure the 6Z(4), 61(5) and other crossings 
in order to establish whether there is a consistent relative 
error between experiment and theory. 

Errors in the level positions can result from deficien- 
cies in either the long-range or the short-range part of the 
model potentials. As discussed above, there are remain- 
ing discrepancies in level positions of up to 0.05 mT for 
s, d, and g levels, and up to 0.5 mT for I levels. Further 
improvements in the potential model are thus needed for 
this important prototype system. This is particularly im- 
portant for predicting the resonances and crossings in the 
80 mT reg ion, where interesting Efimov physics is pre- 
dicted [X 31 ] and even greater sensitivity to model errors is 
expected. A major advantage of the propagator method 
introduced here is that it is inexpensive enough to be 
used to determine model parameters by least-squares fit- 
ting to level energies and locations and strengths of level 
crossings. 

V. CONCLUSIONS 

We have presented a new computational method for 
calculating bound states of molecules such as Cs2. The 
method is based on solving a set of coupled differen- 
tial equations by propagation, without relying on a basis 
set for the interatomic coordinate. This is much more 
efficient than using a radial basis set and allows the 
use of much larger basis sets of spin functions. It also 
eliminates problems with calculating bound states very 
near to dissociation, because the propagation can be ex- 
tended to very large separations at very little expense. 
The new method makes it possible for the first time to 



carry out fully converged calculations on bound states of 
Cs2, including anisotropic couplings due to spin-spin and 
second-order spin-orbit interactions, and to characterize 
avoided crossings between pairs of levels. 

We have compared the results of converged calcula- 
tions using the current best Cs2 model potentials with ex- 
perimental measurements on the near-dissociation states 
of Cs2 in a magnetic field. The model generally per- 
forms well for s, d and g states (with L = 0, 2 and 4), 
though even there there are quantitative discrepancies of 
up to 0.05 mT in the magnetic fields at which levels cross. 
The discrepancies are much larger (0.5 mT) for I states 
(L = 8). The strengths of the avoided crossings also ap- 
pear to be systematically underestimated by the current 
model. These discrepancies should in future allow the 
development of improved models for the potential curves 
and couplings in the Cs2 dimer. Such model improvement 
is both desirable and possible, not only for near-threshold 
levels but also to provide an improved representation of 
more deeply bound states such as those measured by Van- 
haecke et al. [f| . High-quality models are also important 
for proposals to use precision measurements on Cs2 for 
fundamental physics studies [H, H(| . 
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APPENDIX A: MATRIX ELEMENTS 



In the decoupled basis set ([5]), the matrix elements of the isotropic potential operator V C (R) between primitive 
(unsymmetrized) basis functions are 

(sim s ii 1 m il S2m S 2i2m l 2LM L \V c (R)\s 1 m' sl iim' il S2m' s2 i 2 m' i2 L / M' L ) = <^£M z ,M£<Wm< 1 <WK2 

Sl S2 S \ f Si s 2 S 



Yv s (R)(-l) 2si - 2s2+m ° 1+m ° 2+ <i + < 2 (2S + l){ 81 32 a S ) 32 , , . (Al) 

s \ m si rn s2 -m s i-m s2 J \m sl m s2 -m' sl - m s2 ) 

The corresponding matrix elements of the spin-spin operator are 

(sim sl i 1 m il S2m S 2i2^M2LM L \V A (R)\s 1 m' sl iim' il S2m 1 ^ = 5 miim > n 5 mam > a \(R) 

( _ ir+S2 - msl - ms2 -M £ [si(si + 1)(2si + 1)s2(fi2 + 1)(2fi2 + 1)(2X + 1)(2i , + 1}] i/2 2 

V ( L 2 L ' V 1 1 2 )( 31 1 51 )( 32 1 52 ^ (A2) 



M L -q\-qi M' L J \qi q 2 -qi - q%) \-m s i qi m' sl J \-m s2 ?2 m' s2 

where for any individual matrix element the sums over q\ and <j2 collapse because of the selection rules imposed by 
the last two 3-j symbols. The matrix elements of the atomic nuclear spin operators are particularly simple in this 
basis set, 

(sim sl iim il S2m S 2i2m l 2LM L \ii ■ si\sim' sl iim' il s 2 m' s2 i 2 rn' i2 L l M' L ) = 

5LL>fi ms2m ',8 mi2m >{sim, s iiima\ii ■ h Isim'^iim'^} (A3) 
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where 

(siTOgi£imji|?i • si\sim sl i 1 m l i) = m. a m sl ; (A4) 
(sxm sl iima\ii ■ §i\sim sl ± lhma T 1) = [«i(si + 1) - m sl (m sl ± 1)] 1/2 + 1) - TOii(ma =F 1)] 1/2 , (A5) 

and similarly for i2 ■ §2- The matrix elements of L 2 are simply 

(siTO s iiimaS2m S 2«2TOi2iM L |L 2 |siTO' sl iiTO- 1 S2"i(, 2 z 2 m- 2 L'A/^) = 

L(L + 1), (A6) 

and those of the Zccman operator are 

(sim sl iim l iS2m s2 i2'mi2LM L \g e ii- B B s zj + g n fi B B i zj \sim' sl iim[ 1 S2'rn' s2 i2rn' i2 L' M' L ) = 

m sl m sl *^ m il m ii ^ m s2"T-^ 2 ^ rn i% rn i2 

(g e fx B Bm sj + g n fj, B Bmij). (A7) 

All the calculations in the present paper used basis functions symmetrized for exchange of two identical particles 
with si = S2 = s and i\ = i% = i. For m s i = m S 2 or mu = the symmetrized functions are identical to the 
unsymmetrized ones, except that only even L is allowed for bosons and only odd L for fermions. For m s \ ^ m S 2 or 
mu ^ rrii2, the symmetrized functions are 

[\sm s iimiism S 2imi2LM L ) ± (-l) L \sm S 2imi2sm s iimuLM L )] /y/2, (A8) 

with the + sign for bosons and the — sign for fermions. 

The Hamiltonian in the basis set, Eq. (j4|), used in the DVR calculations can be derived from the Hamiltonian in 
the uncoupled basis by performing a unitary transformation, namely, the transformation \oijrrifj) to \sjm S i)\iimn) 
for each of the two atoms (j=l or 2). The transformation depends on the magnetic field strength [371 ]. In practice, 
the eigenvectors for the monomer hj must be evaluated. As rrifj is conserved at most a 2x2 matrix needs to be 
diagonalized. Bose/Fermi symmetrization is ensured by 

[|aiTO/ia2m/2, LM L ) ± (-l) L |a 2 m/ 2 aim/i, LMl)]/V% 

when ai ^ «2 or m/i ^ "i/2- The state with ol\ = 0:2 and m/i = m/2 exists only for even (odd) L for bosonic 
(fermionic) atoms respectively. 
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